packages <- c("CIMseq", "CIMseq.testing", "CIMseq.data", "tidyverse", "ggthemes")
purrr::walk(packages, library, character.only = TRUE)
rm(packages)

##DATA
load('../data/CIMseqData.rda')
load('../data/sObj.rda')
if(!dir.exists('../figures')) dir.create('../figures')
p <- plotUnsupervisedClass(cObjSng, cObjMul, palette('total'))
p

ggsave(
  plot = p,
  filename = '../figures/MGA.enge.classes.pdf',
  device = cairo_pdf,
  height = 180,
  width = 180,
  units = "mm"
)
p <- plotUnsupervisedMarkers(
  cObjSng, cObjMul,
  c("Lgr5", "Ptprc", "Chga", "Dclk1", "Alpi", "Slc26a3", "Atoh1", "Lyz1"),
  pal = RColorBrewer::brewer.pal(8, "Set1")
)
p

ggsave(
  plot = p,
  filename = '../figures/MGA.enge.markers.pdf',
  device = cairo_pdf,
  height = 180,
  width = 180,
  units = "mm"
)

p <- plotUnsupervisedMarkers(cObjSng, cObjMul, "Mki67")
p

ggsave(
  plot = p,
  filename = '../figures/MGA.enge.Mki67.pdf',
  device = cairo_pdf,
  height = 180,
  width = 180,
  units = "mm"
)

# p <- plotUnsupervisedMarkers(cObjSng, cObjMul, "Hoxb13")
# ggsave(
#   plot = p,
#   filename = '../figures/MGA.enge.Hoxb13.pdf',
#   device = cairo_pdf,
#   height = 180,
#   width = 180,
#   units = "mm"
# )

Mouse age

p <- getData(cObjSng, "dim.red") %>%
  matrix_to_tibble("sample") %>%
  inner_join(MGA.Meta, by = "sample") %>%
  mutate(subject_age = as.character(subject_age)) %>%
  mutate(subject_age = parse_factor(subject_age, levels = sort(unique(subject_age)))) %>%
  ggplot() +
  geom_point(aes(V1, V2, colour = subject_age), alpha = 0.7) +
  ggthemes::theme_few() +
  labs(x = "UMAP dim. 1", y = "UMAP dim. 2") +
  scale_colour_manual(values = rev(col40())) +
  guides(colour = guide_legend(title = "Age (days)"))
p

ggsave(
  plot = p,
  filename = '../figures/MGA.enge.QC.age.pdf',
  device = cairo_pdf,
  height = 180,
  width = 180,
  units = "mm"
)

Mouse sex

p <- getData(cObjSng, "dim.red") %>%
  matrix_to_tibble("sample") %>%
  inner_join(MGA.Meta, by = "sample") %>%
  mutate(subject_sex = as.character(subject_sex)) %>%
  mutate(subject_sex = parse_factor(subject_sex, levels = sort(unique(subject_sex)))) %>%
  ggplot() +
  geom_point(aes(V1, V2, colour = subject_sex), alpha = 0.7) +
  ggthemes::theme_few() +
  labs(x = "UMAP dim. 1", y = "UMAP dim. 2") +
  scale_colour_manual(values = rev(col40())) +
  guides(colour = guide_legend(title = "Sex"))

p

ggsave(
  plot = p,
  filename = '../figures/MGA.enge.QC.sex.pdf',
  device = cairo_pdf,
  height = 180,
  width = 180,
  units = "mm"
)

Total counts

cs <- colSums(getData(cObjSng, "counts"))

p <- getData(cObjSng, "dim.red") %>%
  matrix_to_tibble("sample") %>%
  inner_join(tibble(sample = names(cs), counts = cs)) %>%
  ggplot() +
  geom_point(aes(V1, V2, colour = counts), alpha = 0.7) +
  ggthemes::theme_few() +
  labs(x = "UMAP dim. 1", y = "UMAP dim. 2") +
  scale_colour_viridis_c() +
  guides(colour = guide_colourbar(title = "Total counts"))
## Joining, by = "sample"
p

ggsave(
  plot = p,
  filename = '../figures/MGA.enge.QC.counts.pdf',
  device = cairo_pdf,
  height = 180,
  width = 180,
  units = "mm"
)

Total genes

tg <- apply(getData(cObjSng, "counts"), 2, function(c) length(which(c != 0)))

p <- getData(cObjSng, "dim.red") %>%
  matrix_to_tibble("sample") %>%
  inner_join(tibble(sample = names(tg), genes = tg)) %>%
  ggplot() +
  geom_point(aes(V1, V2, colour = genes), alpha = 0.7) +
  ggthemes::theme_few() +
  labs(x = "UMAP dim. 1", y = "UMAP dim. 2") +
  scale_colour_viridis_c() +
  guides(colour = guide_colourbar(title = "Total genes"))
## Joining, by = "sample"
p

ggsave(
  plot = p,
  filename = '../figures/MGA.enge.QC.genes.pdf',
  device = cairo_pdf,
  height = 180,
  width = 180,
  units = "mm"
)

Cell type markers

markers <- c(
  "Lgr5", "Ptprc", "Chga", "Dclk1", "Alpi", "Slc26a3", "Atoh1", "Lyz1", "Mki67",
  "Hoxb13", "Plet1", "Osr2", "Reg4", "Sval1"
)

gene.order <- c(
  "Hoxb13", "Osr2", "Atoh1", "Reg4", "Sval1", "Plet1", "Mki67", "Lgr5", "Alpi", 
  "Slc26a3", "Lyz1", "Dclk1", "Chga", "Ptprc"
)

d <- getData(cObjSng, "counts.cpm")[markers, ] %>%
  t() %>%
  as.data.frame() %>%
  rownames_to_column("Sample") %>%
  as_tibble() %>% 
  gather(gene, cpm, -Sample) %>%
  inner_join(tibble(
    Sample = colnames(getData(cObjSng, "counts")), 
    Classification = getData(cObjSng, "classification")
  ), by = "Sample") %>%
  mutate(Classification = parse_factor(Classification, levels = classOrder.MGA('total'))) %>%
  mutate(gene = parse_factor(gene, levels = gene.order)) %>%
  arrange(Classification, gene) %>%
  mutate(Sample = parse_factor(Sample, levels = unique(Sample))) %>% 
  mutate(type = case_when(
    str_detect(Classification, "^S") ~ "Small intestine",
    str_detect(Classification, "^C") ~ "Colon",
    TRUE ~ "Other"
  )) %>%
  mutate(type = parse_factor(
    type,
    levels = c("Small intestine", "Colon", "Other")
  ))

under <- d %>%
  mutate(idx = 1:nrow(.)) %>% 
  group_by(Classification) %>% 
  summarize(
    min = (min(idx)), max = (max(idx) - 1), median = median(idx)
  ) %>%
  ggplot() +
  geom_segment(
    aes(x = min, xend = max, y = 0, yend = 0, colour = Classification)
  ) +
  geom_text(
    aes(median, y = 0, label = Classification),
    angle = 90, nudge_y = -0.01, hjust = 1, colour = "grey10"
  ) +
  scale_colour_manual(values = palette('total')[classOrder.MGA('total')]) +
  ylim(c(-1, 0)) +
  theme_few() +
  scale_x_continuous(expand = c(0, 0)) +
  theme(
    axis.text = element_blank(),
    axis.title = element_blank(),
    axis.ticks = element_blank(),
    panel.border = element_blank()
  ) +
  guides(colour = FALSE)

over <- d %>% 
  ggplot() + 
  geom_bar(aes(Sample, cpm, fill = Classification), stat = "identity") + 
  facet_grid(gene ~ type, scales = "free", space = "free_x") +
  scale_fill_manual(values = palette('total')[classOrder.MGA('total')]) +
  theme_few() +
  theme(
    axis.text.x = element_blank(),
    axis.title.x = element_blank(),
    axis.ticks.x = element_blank(),
    strip.background.x = element_rect(fill = "white")
  ) +
  labs(y = "CPM") +
  guides(fill = FALSE)

library(patchwork)

p <- over + under + plot_layout(ncol = 1, heights = c(2, 1))
p

ggsave(
  plot = p,
  filename = '../figures/MGA.enge.geneBar.pdf',
  device = cairo_pdf,
  height = 300,
  width = 300,
  units = "mm"
)

Stemness dotplot

markers <- c(
  "Wnt3", "Wnt2b", "Wnt4", "Wnt5a", "Lgr5", "Axin2", "Ascl2", "Egf", "Notum", 
  "Kit", "Dll1", "Dll4", "Reg4", "Plet1"
)
gene.order <- markers

#scale.func <- scale_size
scale.func <- scale_radius
scale.min <- NA
scale.max <- NA
cOrderAlt <- c(
  "SI Goblet", "SI Paneth", "SI Stem", "SI Transit amplifying", 
  "SI Progenitor early", "SI Progenitor late", "SI Enterocytes", 
  "C Goblet proliferating", "C Goblet Plet1 1", "C Goblet Plet1 2", "C Stem 1",
  "C Stem 2", "C Stem 3", "C Transit amplifying", "C Progenitor", 
  "C Colonocytes", "C Goblet Plet1 neg.", "Enteroendocrine", "Tufft", "Blood"
)

p <- getData(cObjSng, "counts.cpm")[markers, ] %>%
  t() %>%
  as.data.frame() %>%
  rownames_to_column("Sample") %>%
  as_tibble() %>% 
  gather(gene, cpm, -Sample) %>%
  inner_join(tibble(
    Sample = colnames(getData(cObjSng, "counts")), 
    Classification = getData(cObjSng, "classification")
  )) %>%
  mutate(type = case_when(
    str_detect(Classification, "^S") ~ "Small intestine",
    str_detect(Classification, "^C") ~ "Colon",
    TRUE ~ "Other"
  )) %>%
  group_by(gene, Classification) %>%
  summarize(mean = mean(cpm), pct = 100 * (length(which(cpm != 0)) / n())) %>%
  mutate(scaled.mean.exp = scale(mean)) %>%
  ungroup() %>%
  mutate(type = case_when(
    str_detect(gene, "^Wnt") ~ "Wnt ligand",
    gene %in% c("Lgr5", "Axin2", "Ascl2") ~ "Wnt target",
    gene %in% c("Egf", "Notum", "Dll1", "Dll4") ~ "Other stemness factors",
    gene %in% c("Plet1", "Reg4", "Kit") ~ "Stem-adjacent\nsecretory cell factors",
    TRUE ~ "NA"
  )) %>%
  mutate(tissue = case_when(
    str_detect(Classification, "^SI") ~ "Small intestine",
    str_detect(Classification, "^C") ~ "Colon",
    TRUE ~ "Other"
  )) %>%
  mutate(tissue = parse_factor(tissue, levels = c("Small intestine", "Colon", "Other"))) %>%
  mutate(Classification = parse_factor(Classification, levels = cOrderAlt)) %>%
  mutate(gene = parse_factor(gene, levels = rev(gene.order))) %>%
  mutate(type = parse_factor(
    type, 
    levels = c(
      "Wnt ligand", "Wnt target", "Other stemness factors", 
      "Stem-adjacent\nsecretory cell factors"
  ))) %>%
  ggplot() +
  geom_point(
    aes(Classification, gene, size = pct, colour = scaled.mean.exp)
  ) +
  facet_grid(type ~ tissue, space = "free", scales = "free") +
  scale_color_gradient(low = "white", high = "darkred") +
  scale.func(range = c(0, 12), limits = c(scale.min, scale.max)) +
  guides(
    size = guide_legend(
      title = "% expressed", title.position = "top", title.hjust = 0.5
    ),
    colour = guide_colorbar(
      title = "Scaled mean expression", title.position = "top", 
      title.hjust = 0.5, barwidth = 10
    )
  ) +
  theme_bw() +
  theme(
    axis.text.x = element_text(angle = 90, hjust = 1), 
    axis.title = element_blank(), 
    legend.position = "top",
    legend.justification = "center"
  )
## Joining, by = "Sample"
p

ggsave(
  plot = p,
  filename = '../figures/MGA.enge.stemness.dotplot.pdf',
  device = cairo_pdf,
  height = 250,
  width = 250,
  units = "mm"
)

Relative frequency of cell types in singlets vs. deconvoluted multiplets

singlets <- c(table(getData(cObjSng, "classification")))
singlets <- singlets / sum(singlets)
deconv <- colSums(adjustFractions(cObjSng, cObjMul, sObj))
deconv <- deconv[match(names(singlets), names(deconv))]
deconv <- deconv / sum(deconv)
if(!identical(names(singlets), names(deconv))) stop("name mismatch")

p <- tibble(
  class = names(singlets),
  singlet.freq = singlets,
  multiplet.freq = deconv
) %>%
  ggplot() +
  geom_point(aes(singlet.freq, multiplet.freq, colour = class), size = 3) +
  scale_colour_manual(values = palette('total')[order(names(palette('total')))]) +
  xlim(min(c(deconv, singlets)), max(c(deconv, singlets))) +
  ylim(min(c(deconv, singlets)), max(c(deconv, singlets))) +
  geom_abline(slope = 1, intercept = 0, lty = 3, colour = "grey") +
  labs(x = "Singlet relative frequency", y = "Multiplet relative frequency") +
  guides(colour = guide_legend(title = "Cell Type")) +
  theme_few()

p

ggsave(
  plot = p,
  filename = '../figures/MGA.enge.sngMulRelFreq.pdf',
  device = cairo_pdf,
  height = 160,
  width = 250,
  units = "mm"
)
plotSwarmCircos(
  sObj, cObjSng, cObjMul, weightCut = 0, classOrder = classOrder.MGA('total'), 
  classColour = palette('total')[classOrder.MGA('total')], h.ratio = 0.85
)
## Joining, by = "class"

Detected duplicates - quadruplicates.
Only ERCC estimated cell number max 4.
Weight cutoff = 10.

# adj <- adjustFractions(cObjSng, cObjMul, sObj, binary = TRUE)
# samples <- rownames(adj)
# rs <- rowSums(adj)
# keep <- rs == 2 | rs == 3 | rs == 4

plotSwarmCircos(
  sObj, cObjSng, cObjMul, weightCut = 10, 
  classOrder = classOrder.MGA('total'), classColour = palette('total')[classOrder.MGA('total')], h.ratio = 0.85,
  alpha = 1e-3
)
## Joining, by = "class"

pdf('../figures/MGA.enge.circos.pdf', width = 12, height = 10, onefile = FALSE)
plotSwarmCircos(
  sObj, cObjSng, cObjMul, weightCut = 10, 
  classOrder = classOrder.MGA('total'), classColour = palette('total')[classOrder.MGA('total')], h.ratio = 0.85,
  alpha = 1e-3
)
## Joining, by = "class"
invisible(dev.off())
sessionInfo()
## R version 3.6.1 (2019-07-05)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Mojave 10.14.6
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] patchwork_0.0.1      ggthemes_4.2.0       forcats_0.4.0       
##  [4] stringr_1.4.0        dplyr_0.8.3          purrr_0.3.2         
##  [7] readr_1.3.1          tidyr_0.8.3          tibble_2.1.3        
## [10] ggplot2_3.2.1        tidyverse_1.2.1      CIMseq.data_0.0.1.4 
## [13] CIMseq.testing_0.0.2 CIMseq_0.3.0.2      
## 
## loaded via a namespace (and not attached):
##  [1] nlme_3.1-141        matrixStats_0.55.0  lubridate_1.7.4    
##  [4] RColorBrewer_1.1-2  httr_1.4.1          gmodels_2.18.1     
##  [7] tools_3.6.1         backports_1.1.4     R6_2.4.0           
## [10] lazyeval_0.2.2      BiocGenerics_0.30.0 colorspace_1.4-1   
## [13] withr_2.1.2         tidyselect_0.2.5    gridExtra_2.3      
## [16] compiler_3.6.1      cli_1.1.0           rvest_0.3.4        
## [19] xml2_1.2.2          labeling_0.3        scales_1.0.0       
## [22] digest_0.6.20       rmarkdown_1.15      pkgconfig_2.0.2    
## [25] htmltools_0.3.6     rlang_0.4.0         GlobalOptions_0.1.0
## [28] readxl_1.3.1        rstudioapi_0.10     shape_1.4.4        
## [31] farver_1.1.0        generics_0.0.2      jsonlite_1.6       
## [34] gtools_3.8.1        magrittr_1.5        Rcpp_1.0.2         
## [37] munsell_0.5.0       S4Vectors_0.22.1    viridis_0.5.1      
## [40] stringi_1.4.3       EngeMetadata_0.1.2  yaml_2.2.0         
## [43] ggraph_2.0.0        MASS_7.3-51.4       plyr_1.8.4         
## [46] Rtsne_0.15          grid_3.6.1          parallel_3.6.1     
## [49] gdata_2.18.0        listenv_0.7.0       ggrepel_0.8.1      
## [52] crayon_1.3.4        lattice_0.20-38     haven_2.1.1        
## [55] graphlayouts_0.5.0  circlize_0.4.8      hms_0.5.1          
## [58] zeallot_0.1.0       knitr_1.24          pillar_1.4.2       
## [61] igraph_1.2.4.1      pso_1.0.3           reshape2_1.4.3     
## [64] future.apply_1.3.0  codetools_0.2-16    stats4_3.6.1       
## [67] glue_1.3.1          evaluate_0.14       modelr_0.1.5       
## [70] vctrs_0.2.0         tweenr_1.0.1        cellranger_1.1.0   
## [73] gtable_0.3.0        RANN_2.6.1          polyclip_1.10-0    
## [76] future_1.14.0       assertthat_0.2.1    xfun_0.9           
## [79] gridBase_0.4-7      ggforce_0.3.1       broom_0.5.2        
## [82] tidygraph_1.1.2     googledrive_1.0.0   viridisLite_0.3.0  
## [85] globals_0.12.4